clear
#delimit ;
set more off;

use pw_coalitions.dta;

xtset warnum;
fp <pceyr>: xtprobit failure maxshift losercinctvc shiftloser lognum loglength <pceyr>, re;

preserve;
set seed 304876538;

drawnorm SN_b1-SN_b9, n(10000) means(e(b)) cov(e(V)) clear;
postutil clear;
postfile mypost prob_hat0 lo0 hi0 prob_hat1 lo1 hi1 diff_hat diff_lo diff_hi
using pwsims3.dta, replace;

noisily display "let's do this shit";

local a=0 ;
while `a' <= 0.33 { ;

{ ;

scalar h_maxshift=0;
scalar h_losercinctvc=0.0569628;
scalar h_number=1.0986123;
scalar h_loglength=5.54;
scalar h_pceyr_1=64;
scalar h_pceyr_2=88;
scalar h_constant=0;

	generate x_betahat0 = SN_b1*h_maxshift
							+ SN_b2*`a'
							+ SN_b3*h_maxshift*`a'
							+ SN_b4*h_number
							+ SN_b5*h_loglength
							+ SN_b6*h_pceyr_1
							+ SN_b7*h_pceyr_2
							+ SN_b8*h_constant;
	generate x_betahat1 = SN_b1*(h_maxshift+0.43)
							+ SN_b2*`a'
							+ SN_b3*(h_maxshift+0.43)*`a'
							+ SN_b4*h_number
							+ SN_b5*h_loglength
							+ SN_b6*h_pceyr_1
							+ SN_b7*h_pceyr_2
							+ SN_b8*h_constant;
	
	gen prob0 =normal(x_betahat0);
	gen prob1 =normal(x_betahat1);
	gen diff=prob1-prob0;
	
	egen probhat0 =mean(prob0);
    egen probhat1=mean(prob1);
    egen diffhat=mean(diff);
    
    tempname prob_hat0 lo0 hi0 prob_hat1 lo1 hi1 diff_hat diff_lo diff_hi;  

    _pctile prob0, p(2.5,97.5) ;
    scalar `lo0' = r(r1);
    scalar `hi0' = r(r2);  
    
    _pctile prob1, p(2.5,97.5);
    scalar `lo1'= r(r1);
    scalar `hi1'= r(r2);  
    
    _pctile diff, p(2.5,97.5);
    scalar `diff_lo'= r(r1);
    scalar `diff_hi'= r(r2);  
   
    scalar `prob_hat0'=probhat0;
    scalar `prob_hat1'=probhat1;
    scalar `diff_hat'=diffhat;
    
    post mypost (`prob_hat0') (`lo0') (`hi0') (`prob_hat1') (`lo1') (`hi1') 
                (`diff_hat') (`diff_lo') (`diff_hi');
   
    };      
    
    drop x_betahat0 x_betahat1 prob0 prob1 diff probhat0 probhat1 diffhat;
    
    local a=`a'+ 0.001;
    
    display "." _c;
    
} ;

display "";

postclose mypost;

restore;

clear;

use pwsims3.dta;

gen MV = _n-1;
replace MV=MV/1000;

twoway (line diff_hat MV if MV<=0.3, sort) (line diff_lo MV if MV<=0.3, sort lpattern(dash)) (line diff_hi MV if MV<=0.3, sort lpattern(dash)), ytitle("Change in Pr(War), Min to Max Shift") yline(0, lcolor(gray)) xtitle("Losing Side's Total Military Capabilities") legend(off);
graph save Graph figure5b.gph, replace;

save pwsims3.dta, replace;
